program define levels2_byhc
   version 15.1	
	* future: extend to allow label1 label2 and legend
	* revised 1/30/2021
	syntax varlist [if] [pweight], GP1(varname) GP2(varname) [OPTS(str) ALLFIVE BYOPTS(str)]
	marksample touse
	
	if "`weight'" ~= "" {
	   local wt [`weight'`exp']
   }	
	
   forvalues h = 0/1 {
	   di _n(2) "** Humans causing `h' **"
      levels2 `varlist' if `touse' & hc2==`h' `wt', gp1(`gp1') gp2(`gp2') ///
		   opts(`opts') `allfive' nodraw
		tempname lvl1`h' lvl2`h'
      mat `lvl1`h'' = r(lvl1)
      mat `lvl2`h'' = r(lvl2)
   }

	* combine plots horizontally
	local gopts scheme(s1mono) plotregion(style(none)) mlabel ///
	   format(%3.0f) xlabel(0(20)100)
	coefplot (matrix(`lvl11'[1,]), ci("`lvl11'[5,] `lvl11'[6,]") msymbol(O) mlabpos(2))  ///
		(matrix(`lvl21'[1,]), ci("`lvl21'[5,] `lvl21'[6,]") msymbol(Oh) mlabpos(10) /// 
		mcolor(gs8) mlabcolor(gs8)), bylabel("Humans Causing") ///		
		||  ///
	   (matrix(`lvl10'[1,]), ci("`lvl10'[5,] `lvl10'[6,]") msymbol(O) mlabpos(2))  ///
		(matrix(`lvl20'[1,]), ci("`lvl20'[5,] `lvl20'[6,]") msymbol(Oh) mlabpos(10) /// 
		mcolor(gs8) mlabcolor(gs8)), bylabel("Not Causing or DK") ///
		coeflabel(none = "No Action" low = "Cut 5%" high = "Cut 25%") /// only applies if not all 5
		`gopts' `opts' ///
		byopts(rows(1) legend(off) `byopts') subtitle(, fcolor(white))
	
end
